#### MAKE LAGS BASED ON FACTORS

lag <- function(x, s, factor=NULL){
  lagtemp <- function(x,s) {
    if(length(x) > s)  M = c(rep(NA, length=s), rev(rev(x)[-(1:s)]))
    if(length(x) <= s) M = rep(NA, length=length(x))
    M
  }
  lagtemp2 <- function(x, s, factor){
    if (is.null(factor)) return(lagtemp(x, s))
    if (!is.null(factor)) return(unlist(tapply(x, factor(factor), lagtemp, s=s)))
  }
  sapply(1:s, lagtemp2, x=x, factor=factor)
}


##### Combine results from regressions

##### Combine results from regressions

outtex <- function(models, below = TRUE, ci = NULL, omit = NULL, only = NULL, single.var = FALSE, caption = NULL, label = NULL, var.labels = NULL, model.names = NULL, col.names = "", table.placement = 't!', space = "50pt", robust = FALSE, size = "normalsize", add.line = list(pos = NULL, command = NULL), file = ""){
  
  outr <- function(model, below, ci, omit, only, single.var){
    
    output <- data.frame(coefficients(summary(model)))
    
    if (!is.null(omit)) vars <- which(!grepl(paste(omit, collapse="|"), rownames(output)))
    
    if (!is.null(only)) vars <- which(grepl(paste(omit, collapse="|"), rownames(output)))
    
    if(robust) output <- coeftest(model,  vcov=function(x) vcovHC(x, type="HC1"))
    
    output  <- output[vars, ]
    varnames <- rownames(output)
    
    if(class(model)[1]%in%c("lmerMod", "glmerMod")) output$p <- 2*pnorm(abs(output[, 3]), lower.tail = FALSE)
    
    r <- output[, 1:2]
    p <- output[, 4]
    p <- ifelse(p < 0.001, "^{\\dagger}", ifelse(p < 0.01, "^{**}", ifelse(p < 0.05, "^{*}", "")))
    se <- paste("(", roundr(r[,2], 2), ")", sep="")
    
    if (!is.null(ci)) {
      ci$method <- ifelse(!is.list(ci), "profile", ci$method)
      ci <- roundr(confint(model, par = varnames, method = ci$method), 2)
      se <- paste("(", ci[,1], ", ", ci[,2], ")" , sep = "")
      p <- ""
    }
    
    m <- data.frame(names = varnames, b = paste(roundr(r[,1]), p, sep=""), se = se, stringsAsFactors = FALSE)
    
    if (class(model)[1]%in%c("lm", "lmrob")) {
      r2 <- roundr(summary(model)$adj.r.squared, 2)
      names(r2) <- "Adj. $R^2$"
    }
    
    if (class(model)[1]%in%c("plm")) {
      r2 <- roundr(summary(model)$r.squared[2])
      names(r2) <- "Adj. $R^2$"
    }
    
    if(class(model)[1]%in%c("lmerMod")) {
      r2 <- roundr(cor(predict(model), model.frame(model)[,1])^2, 2)
      names(r2) <- "Adj. $R^2$"
    }
    
    if (class(model)[1]%in%c("glm", "glmerMod")){
      r2 <- prettyNum(roundr(logLik(model), 0), big.mark=",", preserve.width="none")
      names(r2) <- "Log-likelihood"
    }
    
    if (class(model)[1]=="rlm") {
      r2 <- roundr(cor(predict(model), model.frame(model)[,1])^2,2)
      names(r2) <- "Adj. $R^2$"
    }
    
    n <- prettyNum(nrow(model.matrix(model)), big.mark=",", preserve.width="none")      
    
    if(below==FALSE) { 
      if(single.var == TRUE)  m <- data.frame(m, r2 = r2, n = n)
      if(single.var == FALSE) {
        fit <- c(names(r2), paste("\\multicolumn{1}{r}{", r2, "}", sep = ""), "")
        n <- c("Observations", paste("\\multicolumn{1}{r}{", n, "}", sep = ""), "")
        m <- rbind(m, fit, n)
        m <- data.frame(m)
      }}
    
    if(below==TRUE)  {
      
      r2 <- paste("\\multicolumn{1}{c}{", r2, "}", sep = "")
      n <- paste("\\multicolumn{1}{c}{", n, "}", sep = "")  
      est <- c(t(m[,-1]))
      est <- as.matrix(c(est, r2, n))
      names <- c(t(cbind(m$names, paste(m$names, "se", sep = "-"))))
      names <- c(names, "R2", "N")
      names[1] <- "Intercept"
      if (class(model)[1]=="glm") names[which(names=="R2")] <- "Log-Likelihood"
      
      m <- data.frame(names = names, est = est, stringsAsFactors = FALSE)
    }
    
    list(m = m, varnames = varnames)
  }
  
  rr <- lapply(models, outr, ci = ci, below = below, omit = omit, only = only, single.var = single.var)
  varnames <- unique(unlist(lapply(rr, function(x) x$varnames)))
  rr <- lapply(rr, function(x) x$m)
  m <- suppressWarnings(Reduce(function(...) merge(..., by = "names", all=TRUE), rr))
  rownames(m) <- m$names
  
  sort.names <- c(varnames, setdiff(Reduce(union, lapply(rr, function(x) x$names)), varnames))
  
  #sort.names <- c(sort.names, rr[[1]]$names[-(1:(length(x$names)-2))])
  
  if (!below) {
    m <- m[sort.names, ]
    if(!is.null(var.labels)) {
      m[varnames, "names"] <- var.labels
    }
    if(is.null(model.names)) model.names <- paste(col.names, 1:((ncol(m)-1)/2))
    return(
      print(xtable(m, 
                   caption = caption,
                   label = label,
                   align = c("l@{\\hskip -10pt}", "l", paste(rep("d", ncol(m)-2), rep(c("@{\\hskip -7pt}", paste("@{\\hskip ", space, "  }", sep ="")), (ncol(m)-1)/3), sep = ""), "d@{\\hskip 20pt}")), 
            sanitize.text.function = identity,
            size = size,
            file = file,
            hline.after = NULL,
            include.rownames = FALSE, 
            include.colnames = FALSE,
            table.placement = table.placement,
            caption.placement = 'bottom',
            add.to.row = list(pos = as.list(c(add.line$pos, -1, 0, 0, nrow(m)-2, nrow(m), nrow(m))), 
                              command = c(add.line$command, "\\toprule\n", paste("&", paste(paste("\\multicolumn{2}{c}{", model.names, "}", sep = ""), collapse = "&"), "\\\\\n",collapse=""), "\\midrule\n", "\\midrule\n", "\\bottomrule\n", ifelse(!is.null(ci), paste("\\multicolumn{", ncol(m), "}{l}{{\\footnotesize Note: 95 percent confidence intervals in parentheses.\n}}", sep = "") , paste("\\multicolumn{", ncol(m), "}{l}{{\\footnotesize Note:", ifelse(robust, " Robust standard ", " Standard "),"errors in parentheses; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{\\dagger}$p$<$0.001.\n}}", sep = ""))))
      ))
  }
  
  if (below == TRUE){
    
    colnames(m) <- c("names", paste("Model", 1:length(rr), sep = " "))
    
    # Put the stats at the end of the table
    t <- which(sort.names%in%c("R2", "N"))
    sort.names <- c(sort.names[-t], "R2", "N")
    
    m <- m[sort.names, ]
    
    t <- which(grepl("-se", m$names))
    m$names[t] <- ""
    
    m$names[-t]
    
    # Variable names/labels
    
  }
  m[rownames(m)!="NA", ]
  
}




# Hausman test for lmer objects
phtest_glmer <- function (Mod1, Mod2, ...)  {  
  coef.1 <- coefficients(summary(Mod1))[,1]
  coef.2 <- coefficients(summary(Mod2))[,1]  
  vcov.1 <- vcov(Mod1)
  vcov.2 <- vcov(Mod2)
  names.1 <- names(coef.1)
  names.2 <- names(coef.2)
  coef.h <- names.2[names.2 %in% names.1]
  dbeta <- coef.1[coef.h] - coef.2[coef.h]
  df <- length(dbeta)
  dvcov <- vcov.2[coef.h, coef.h] - vcov.1[coef.h, coef.h]
  stat <- abs(t(dbeta)%*%as.matrix(solve(dvcov)) %*% dbeta)  ## added as.matrix()
  pval <- pchisq(stat, df = df, lower.tail = FALSE)
  names(stat) <- "chisq"
  parameter <- df
  names(parameter) <- "df"
  alternative <- "one model is inconsistent"
  res <- list(statistic = stat, p.value = pval, parameter = parameter, 
              method = "Hausman Test",  alternative = alternative)
  class(res) <- "htest"
  res
}


# Nice rounding of numbers for output
roundr <- function(x, n = 2) {
  f <- function(x, n) sprintf(paste("%.", n, "f", sep=""), round(x, n))
  if (is.null(dim(x))) out <- f(x, n)
  if (!is.null(dim(x))) out <- apply(x, 2, f, n = n)
  out
}

#### CALCULATE SURVIVAL TIMES

surv <- function(x){

    if(length(x)==0){return("No observations")}
  
  t=which(x==1)
  tau=1:length(x)
  
  if(length(t)==0)return(tau)
  if(length(t)==1 & x[length(x)]==1)return(tau)
  if(length(t)==1 & x[length(x)]==0)return(c(1:t,1:(length(x)-t)))
  
  if (x[1]==0 & x[length(x)]==0) t=c(1,t,length(x))
  if (x[1]==1 & x[length(x)]==0) t=c(t,length(x))
  if (x[1]==0 & x[length(x)]==1) t=c(1,t)
  
  tau=1:t[2]
  for (i in 3:length(t)){
    tau=c(tau,1:(t[i]-t[i-1]))
  }
  return(tau)
}


surv.until <- function(x, jump=FALSE){
  
  t=which(x==1)
  L=length(x)
  l=length(t)
  
  if(l==0) out = paste(L:1, "+", sep="")
  
  if(l > 0) {
    
    if (jump==TRUE){
      t = t-1
      T = t - c(0, t[-l])
      T = T[T!=0]
      out = c(unlist(sapply(T, function(x) rev(rep(1:x)))))
      out = c(out, paste((L-length(out)):1,"+", sep=""))      
    }
    
    if(jump==FALSE){
      T = t - c(0, t[-l])
      out = c(unlist(sapply(T, function(x) rev(rep(1:x)-1))))
      
      if(x[L]==0) out = c(out, paste((L-length(out)):1,"+", sep=""))
    }
  }
  
  out
  
}


indicate <- function(x){ #create factor to indicate different periods  
  
  t=which(x==1)
  L=length(x)
  l=length(t)
  
  if(l==0) out = rep(1, L)
  
  if (l > 0){
    T = t - c(0, t[-l])
    tau=1:l
    out = unlist(sapply(tau, function(x) rep(tau[x],each=T[x])))
    
    if(x[L]==0){
      out = c(out, rep(l+1, L-length(out)))
    }
  }
  
  out
  
}


